Project 3

1)Reading grey photos and determining the significance level

Firstly we turned the images into greyscale format to use our glcm function. Then, we read the images with raster function. To calculate the control limits, we decided to take significance level 0.003.

library(raster)
## Warning: package 'raster' was built under R version 3.6.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.6.2
library(glcm)
## Warning: package 'glcm' was built under R version 3.6.2
photo_1 <-raster("C:/Users/user/Desktop/Images/Fabric1g.jpg")
photo_2 <-raster("C:/Users/user/Desktop/Images/Fabric2g.jpg")
photo_3 <-raster("C:/Users/user/Desktop/Images/Fabric3g.jpg")
photo_4 <-raster("C:/Users/user/Desktop/Images/Fabric4g.jpg")
photo_5 <-raster("C:/Users/user/Desktop/Images/Fabric5g.jpg")
photo_6 <-raster("C:/Users/user/Desktop/Images/Fabric6g.jpg")
photo_7 <-raster("C:/Users/user/Desktop/Images/Fabric7g.jpg")
photo_8 <-raster("C:/Users/user/Desktop/Images/Fabric8g.jpg")
photo_9 <-raster("C:/Users/user/Desktop/Images/Fabric9g.jpg")
photo_10 <-raster("C:/Users/user/Desktop/Images/Fabric10g.jpg")
photo_11 <-raster("C:/Users/user/Desktop/Images/Fabric11g.jpg")
photo_12 <-raster("C:/Users/user/Desktop/Images/Fabric12g.jpg")
photo_13 <-raster("C:/Users/user/Desktop/Images/Fabric13g.jpg")
photo_14 <-raster("C:/Users/user/Desktop/Images/Fabric14g.jpg")
photo_15 <-raster("C:/Users/user/Desktop/Images/Fabric15g.jpg")
photo_16 <-raster("C:/Users/user/Desktop/Images/Fabric16g.jpg")
photo_17 <-raster("C:/Users/user/Desktop/Images/Fabric17g.jpg")
photo_18 <-raster("C:/Users/user/Desktop/Images/Fabric18g.jpg")
photo_19 <-raster("C:/Users/user/Desktop/Images/Fabric19g.jpg")
photo_20 <-raster("C:/Users/user/Desktop/Images/Fabric20g.jpg")

significance_level=0.003

2) Applying GLCM method and finding Control Limits

In the GLCM function, we selected window size as 11x11 because too small window size make the data much bigger and it looks like monitoring pixel by pixel. On the other hand, larger window sizes cause data loss. We can’t see outliers clearly.

Additionally we selected our grey scale level as 32.If all 256 x 256 cells were used, there would be many cells filled with 0’s (because that combination of grey levels simply does not occur on the image). The GLCM approximates the joint probability distribution of two pixels. Having many 0’s in cells makes this a very bad approximation.

After finding our result, we converted it as a matrix and then calculate the control limits for GLCM mean. In addition, only for the 1st image, we marked the outliers with black.

result1 =glcm(photo_1, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result1$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.559173417068062" "upper bound=0.763181991433118"
m[is.na(m)] = 0
for(i in 6:506){
  for(j in 6:506){
    if(m[i,j]>upper_bound){
      m[i,j]=0
    }
    
  }
  
}

plot(0.5:512.5,0.5:512.5,type="n", main="Image with Defects", xlab = "Horizontal", ylab = "Vertical")
rasterImage(m, 0.5,0.5, 512.5, 512.5,interpolate=FALSE)

plot(result1)

3) Monitoring all images

We applied the GLCM function to all other images to find the outliers.

result2 =glcm(photo_2, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result2$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.450554877988729" "upper bound=0.647530075554598"
plot(result2)

result3 =glcm(photo_3, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result3$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.65619866488857"  "upper bound=0.754297021940252"
plot(result3)

result4 =glcm(photo_4, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result4$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.594660273108069" "upper bound=0.740258605044717"
plot(result4)

result5 =glcm(photo_5, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result5$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.489308479310274" "upper bound=0.595811583766053"
plot(result5)

result6 =glcm(photo_6, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result6$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.261825025335383" "upper bound=0.426047204378742"
plot(result6)

result7 =glcm(photo_7, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result7$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.521469675965352" "upper bound=0.719432411388866"
plot(result7)

result8 =glcm(photo_8, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result8$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.52129296499106"  "upper bound=0.686871796137509"
plot(result8)

result9 =glcm(photo_9, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result9$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.286355572218554" "upper bound=0.378289809538593"
plot(result9)

result10 =glcm(photo_10, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result10$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.369175549275657" "upper bound=0.446106953126564"
plot(result10)

result11 =glcm(photo_11, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result11$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.590756690965596" "upper bound=0.696600756072147"
plot(result11)

result12 =glcm(photo_12, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result12$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.0981017742757438" "upper bound=0.316715628625067"
plot(result12)

result13 =glcm(photo_13, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result13$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.409457176645266" "upper bound=0.581787857305428"
plot(result13)

result14 =glcm(photo_14, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result14$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.511820361872029" "upper bound=0.576420340865059"
plot(result14)

result15 =glcm(photo_15, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result15$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.392134475536488" "upper bound=0.556101419151409"
plot(result15)

result16 =glcm(photo_16, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result16$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.665331156025305" "upper bound=0.890730102397932"
plot(result16)

result17 =glcm(photo_17, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result17$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.23519726928005"  "upper bound=0.415081538806257"
plot(result17)

result18 =glcm(photo_18, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result18$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.472243106441"    "upper bound=0.621767635247112"
plot(result18)

result19 =glcm(photo_19, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result19$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.372698159690963" "upper bound=0.591117637176047"
plot(result19)

result20 =glcm(photo_20, n_grey=32 ,
              window = c(11,11), 
              shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
              statistics = c("mean", "variance", "correlation", "homogeneity", "contrast","dissimilarity")
)

m=as.matrix(result20$glcm_mean)
mean=mean(m,na.rm=TRUE)
sd=sd(m,na.rm=TRUE)
upper_bound = mean + qnorm((1-significance_level/2))*sd 
lower_bound = mean + qnorm(significance_level/2)*sd
print(c(paste0("lower bound=",lower_bound),paste0("upper bound=",upper_bound)))
## [1] "lower bound=0.456961806090594" "upper bound=0.577583674229758"
plot(result20)

```